// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#ifndef PARAM_INTERPRET_H
#define PARAM_INTERPRET_H

// standard headers
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <iomanip>
#include <cmath>
#include <cstdarg>
#include <stdexcept>
#include <utility>
#include <tuple>

#include <Eigen/Eigen>
#include <cmdioutil.h>
#include "xyz.h"


// this went too complex... but it is usefull --> normal grammar for vector expression needed
Eigen::Vector3d interpret_coord(const xyz & s,const std::string & e)
{       
        Eigen::Vector3d v(0.0,0.0,0.0);
        util::triple_expression<double,double,double,'/'> ex;
        util::pair_expression<double,std::string,'*'> exvar; 
        util::pair_expression<int,int,'%'> exii; 
        util::triple_expression<std::string,int,int,'%'> exsii; 
        util::pair_expression<std::string,std::string,'+'> exsum; 
        int index; 
        if (e == "a1") v = s.lattice[0];
        else if (e == "a2") v = s.lattice[1];
        else if (e == "a3") v = s.lattice[2];
        else if (e == "na1") v = s.lattice[0].normalized();
        else if (e == "na2") v = s.lattice[1].normalized();
        else if (e == "na3") v = s.lattice[2].normalized();
        else if (e == "-a1") v = -s.lattice[0];
        else if (e == "-a2") v = -s.lattice[1];
        else if (e == "-a3") v = -s.lattice[2];
        else if (e == "-na1") v = -s.lattice[0].normalized();
        else if (e == "-na2") v = -s.lattice[1].normalized();
        else if (e == "-na3") v = -s.lattice[2].normalized();
        else if (e == "x") v = Eigen::Vector3d(1.,0.,0.);
        else if (e == "y") v = Eigen::Vector3d(0.,1.,0.);
        else if (e == "z") v = Eigen::Vector3d(0.,0.,1.);
        else if (e == "-x") v = Eigen::Vector3d(-1.,0.,0.);
        else if (e == "-y") v = Eigen::Vector3d(0.,-1.,0.);
        else if (e == "-z") v = Eigen::Vector3d(0.,0.,-1.);
        else if(util::to(e,exsum))
        {
                std::string es1 = std::get<0>(exsum); 
                std::string es2 = std::get<1>(exsum); 
                Eigen::Vector3d v1 = interpret_coord(s,es1);
                Eigen::Vector3d v2 = interpret_coord(s,es2);
                v = v1+v2;
        }
        else if(util::to(e,exvar))
        {
                std::string evar = std::get<1>(exvar); 
                if (evar == "a1") v = s.lattice[0];
                else if (evar == "a2") v = s.lattice[1];
                else if (evar == "a3") v = s.lattice[2];
                else if (evar == "na1") v = s.lattice[0].normalized();
                else if (evar == "na2") v = s.lattice[1].normalized();
                else if (evar == "na3") v = s.lattice[2].normalized();
                else if (evar == "-a1") v = -s.lattice[0];
                else if (evar == "-a2") v = -s.lattice[1];
                else if (evar == "-a3") v = -s.lattice[2];
                else if (evar == "-na1") v = -s.lattice[0].normalized();
                else if (evar == "-na2") v = -s.lattice[1].normalized();
                else if (evar == "-na3") v = -s.lattice[2].normalized();
                else if(util::to(evar,exsii)) 
                {
                        std::string evar2 = std::get<0>(exsii);
                        int i1 = convert_index(s.n,std::get<1>(exsii));  
                        int i2 = convert_index(s.n,std::get<2>(exsii));
                        Eigen::Vector3d v1 = s.coords[i1];
                        Eigen::Vector3d v2 = s.coords[i2];
                
                        if (evar2 == "x")      v << v2.x()-v1.x(),0,0;
                        else if (evar2 == "y") v << 0,v2.y()-v1.y(),0;
                        else if (evar2 == "z") v << 0,0,v2.z()-v1.z();
                        else throw std::runtime_error("Unknown displacement constant in const expression: "+e);
                }
                else throw std::runtime_error("Unknown coord constant in coord expression: "+e);
                v *= std::get<0>(exvar);
        }
        else if(util::to(e,exii))
        {
                int i1 = convert_index(s.n,std::get<0>(exii));
                int i2 = convert_index(s.n,std::get<1>(exii));
                Eigen::Vector3d v1 = s.coords[i1];
                Eigen::Vector3d v2 = s.coords[i2];
                v = v2-v1;
        }
        else if(util::to(e,exsii)) 
        {
                std::string evar = std::get<0>(exsii);
                int i1 = convert_index(s.n,std::get<1>(exsii));  
                int i2 = convert_index(s.n,std::get<2>(exsii));
                Eigen::Vector3d v1 = s.coords[i1];
                Eigen::Vector3d v2 = s.coords[i2];
                
                if (evar == "x")      v << v2.x()-v1.x(),0,0;
                else if (evar == "y") v << 0,v2.y()-v1.y(),0;
                else if (evar == "z") v << 0,0,v2.z()-v1.z();
                else throw std::runtime_error("Unknown displacement constant in const expression: "+e);
        }
        else if(util::to(e,ex)) triple2array(ex,v);
        else if(util::to(e,index)) v = s.coords[convert_index(s.n,index)];
        else throw std::runtime_error("Wrong coord expression: "+e);
        return v;
}

double interpret_angle(const xyz & s,const std::string & e)
{
        util::triple_expression<int,int,int,'%'> exiii;
        double angle = 0.0;                                                                                                                                                                                                                                                                                                        
        if(util::to(e,exiii))
        {
                int i1 = convert_index(s.n,std::get<0>(exiii));
                int i2 = convert_index(s.n,std::get<1>(exiii));
                int i3 = convert_index(s.n,std::get<2>(exiii));
                Eigen::Vector3d v1 = s.coords[i1];
                Eigen::Vector3d v2 = s.coords[i2];
                Eigen::Vector3d v3 = s.coords[i3];
                angle = acos((v1-v2).dot(v3-v2)/(v1-v2).norm()/(v3-v2).norm());
        } 
        else if(util::to(e,angle)) angle=angle/180*3.141592653589793;
        else throw std::runtime_error("Wrong angle expression: "+e);
        return angle; 
}

double interpret_disp(const xyz & s,const std::string & e)
{
        util::pair_expression<double,std::string,'*'> exvar; 
        util::triple_expression<std::string,int,int,'%'> exsii; 
        double disp;                                                                                                                                                                                                                                                                                                        
        if (e == "a1") disp = s.lattice[0].norm();
        else if (e == "a2") disp = s.lattice[1].norm();
        else if (e == "a3") disp = s.lattice[2].norm();
        else if (e == "-a1") disp = -s.lattice[0].norm();
        else if (e == "-a2") disp = -s.lattice[1].norm();
        else if (e == "-a3") disp = -s.lattice[2].norm();
        else if(util::to(e,exvar))
        {
                std::string evar = std::get<1>(exvar);
                if (evar == "a1") disp = s.lattice[0].norm();
                else if (evar == "a2") disp = s.lattice[1].norm();
                else if (evar == "a3") disp = s.lattice[2].norm();
                else if (evar == "-a1") disp = -s.lattice[0].norm();
                else if (evar == "-a2") disp = -s.lattice[1].norm();
                else if (evar == "-a3") disp = -s.lattice[2].norm();
                else if(util::to(e,exsii)) 
                {
                        std::string evar2 = std::get<0>(exsii);
                        int i1 = convert_index(s.n,std::get<1>(exsii));
                        int i2 = convert_index(s.n,std::get<2>(exsii));
                        Eigen::Vector3d v1 = s.coords[i1];
                        Eigen::Vector3d v2 = s.coords[i2];
                        
                        if (evar2 == "x")     disp = v2.x()-v1.x();
                        else if (evar2 == "y") disp = v2.y()-v1.y();
                        else if (evar2 == "z") disp = v2.z()-v1.z();
                        else if (evar2 == "d") disp = (v2-v1).norm();
                        else throw std::runtime_error("Unknown displacement constant in displacement expression: "+e);
                }
                else throw std::runtime_error("Unknown coord constant in coord expression: "+e);
                disp *= std::get<0>(exvar);
        }
        else if(util::to(e,exsii)) 
        {
                std::string evar = std::get<0>(exsii);
                int i1 = convert_index(s.n,std::get<1>(exsii));
                int i2 = convert_index(s.n,std::get<2>(exsii));
                Eigen::Vector3d v1 = s.coords[i1];
                Eigen::Vector3d v2 = s.coords[i2];
                
                if (evar == "x")     disp = v2.x()-v1.x();
                else if (evar == "y") disp = v2.y()-v1.y();
                else if (evar == "z") disp = v2.z()-v1.z();
                else if (evar == "d") disp = (v2-v1).norm();
                else throw std::runtime_error("Unknown displacement constant in displacement expression: "+e);
        }
        else if(util::to(e,disp)) ;
        else throw std::runtime_error("Wrong displacement expression: "+e);

        return disp; 
}

std::string interpret_expression_help()
{
        std::string help ="";
        help += "coord expr @ : [<factor>*](a1|a2|a3|-a1|-a2|-a3|na1|na2|na3|-na1|-na2|-na3) | \n";
        help += "               x|y|z|-x|-y|-z | <index>\%<index> | <index> | <x>/<y>/<z> |\n";
        help += "               [<factor>*](x|y|z)\%<index>\%<index>\n";
        help += "displ expr & : [<factor>*](a1|a2|a3|-a1|-a2|-a3) |\n";
        help += "               [<factor>*](x|y|z|d)\%<index>\%<index> | <displacement>\n";
        help += "angle expr $ : <index>\%<index>\%<index> | <degree>\n";
        help += "index expr i : <index>\n";
        return help;
}


// set expression needed (set of atoms)
xyz& select_by_expression(xyz & s,const std::string & e)
{
        util::list_expression<int> l;
        
        for(int i = 0; i < s.n; i++)
        if( e == s.labels[i])
        s.selected[i] = true;
        
        if(util::to(e,l)) select(s,l);

        return s; 
}


#endif // PARAM_INTERPRET_H
